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Formation and competition of associations are studied in a six-species ecological model where each species 
has two predators and two prey. Each site of a square lattice is occupied by an individual belonging to one of 
the six species. The evolution of the spatial distribution of species is governed by iterated invasions between the 
neighboring predator-prey pairs with species specific rates and by site exchange between the neutral pairs with 
a probability X. This dynamical rule yields the formation of five associations composed of two or three species 
with proper spatiotemporal patterns. For large X a cyclic dominance can occur between the three two-species 
associations whereas one of the two three-species associations prevails in the whole system for low values of 
X in the final state. Within an intermediate range of X all the five associations coexist due to the fact that 
cyclic invasions between the two-species associations reduce their resistance temporarily against the invasion of 
three-species associations. 



PACS numbers: 87.23.Cc, 89.75.Fb, 05.50.+q 

I. INTRODUCTION 

Many real systems consist of small different objects whose 
organization into large spatial associations (communities) is 
the result of some evolutionary rules controlling the system's 
behavior at the microscopic level [[ilSlllHlSiifi]. At a larger 
spatial scale the mentioned associations can be considered as 
objects forming larger (super) associations and the repetition 
of this process can even yield a hierarchy of associations. Now 
some general and elementary features of this complex process 
are revealed by a toy model exhibiting several ways how the 
associations coexist. 

The spatial predator-prey models with many species proved 
to be an appropriate model to study the formation and compe- 
tition of associations |7, 8, 9]. In these models the associations 
are composed of a portion of all the species and are character- 
ized by a spatio-temporal pattern. In fact, the associations 
are possible solutions and some of these solutions can be ob- 
served as a final state when the numerical simulations are per- 
formed on small systems. As the solutions of any subsystem 
(where several species are missing) are also solutions for the 
whole system therefore the number of solutions (possible as- 
sociations) increases exponentially with the number of species 
(excepting for some particular food webs). In some cases, in 
spite of the large number of possible solutions, the evolution- 
ary process selects one of the possible solutions characteriz- 
ing the final stationary state even for an infinitely large system 
size. In other cases, equivalent associations compete for ter- 
ritories through a domain growing process, as it happensj^or 
the g-state Potts model below the critical temperature jlOtl . 
and finally one of the associations will prevail in the whole 
(finite) system. Within a wide range of parameters, however, 
the domain growing process is stopped and one can observe a 
self-organizing domain structure (sustaining all species alive) 
where large domains of associations can be clearly identified. 
The self-organizing pattern can be maintained by cyclic dom- 
inance between the associations or by other dynamical phe- 
nomena (sometimes resembling the death and rebirth of the 
Phoenix bird) where different length- and time-scales emerge 



(for examples see the Refs. |H,|^). 

Now we describe another mechanism yielding a self- 
organizing pattern with five associations representing two ba- 
sically different classes of the defensive alliances which can 
be considered as privileged associations. This effect is ob- 
served in a six-species predator-prey model which is a sim- 
plified combination of two previously investigated models 

MM. 



II. THE MODEL 

We consider a six-species predator-prey model where each 
site i of a square lattice is occupied by an individual belonging 
to one of the six species. The species distribution is charac- 
terized by the set of site variables (si = 0, . . . , 5) referring to 
the label of species at the given site i. The predator-prey re- 
lations are defined by a food web indicating that each species 
has two predators and two prey. For the present model we 
distinguish two invasion rates, a and 7 (0 < a, 7 < 1), as 
demonstrated in Fig. [T] The different values of a and 7 pa- 
rameters describe the cases when the strengths of dominance 
within a cyclic alliance and between the members of different 
alliances are unequal. 

The evolution of species distribution is controlled by re- 
peating the following elementary steps. First, two neighbor- 
ing sites (i and j) are chosen at random. If Si is the predator 
of Sj then the site j will be occupied by an offspring of the 
species Si (in short, Sj Si) with a probability given by the 
coiTesponding invasion rate (a or 7). Evidently, for the oppo- 
site predator-prey relation Si will be transformed to the state 
Sj (si Sj) with the corresponding probability. If Si and Sj 
are a neutral pair (e.g., Si = and Sj = 3) then they exchange 
their sites [(si, Sj) — > {sj, Si)] with a probability X charac- 
terizing the strength of mixing. Finally, nothing happens if 

Si ~ Sj. 

The system is started from a random initial state where each 
species is present with the same probability. After many rep- 
etition of the above elementary steps the system evolves into 
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FIG. 1: Food web for the present six-species predator-prey model. 
Arrows point from a predator towards its prey with heterogeneous 
invasion rates specified along the edges. 



a state that can be characterized by the average densities ps 
is ~ 0, . . . , 5) satisfying the condition X]s=o P<> ~ ^- 
many cases the quantification of the nearest-neighbor pair cor- 
relations is necessary to give an adequate description about the 
spatial distribution of species. Therefore, we can introduce 
four types of pair configuration probabilities for the present 
model: pid denotes the probability of finding identical species 
on two neighboring sites; p„ is the probability of finding a 
neutral pair (e.g., species and 3); Pa and p^ are the sum of 
those predator-prey pair probabilities where the invasion rates 
are a and 7, respectively. These quantities are also satisfying 
a normalization condition, i.e., pid + Pn + + p-y = 1. 

The above system was investigated by Monte Carlo (MC) 
simulations on a square lattice of size N = L x L under peri- 
odic boundary conditions and the linear size L is varied from 
400 to 4000. The MC simulations were performed systemat- 
ically for a fixed value of the highest invasion rate (e.g., for 
7 = 1) while the other invasion rate and X are varied grad- 
ually. The stationary states were characterized by the above 
mentioned order parameters averaged over a sampling ts af- 
ter a suitable thermalization time tth- To observe the actual 
spatio-temporal pattern at a specified values of a, 7, and X, 
the parameters L, ts, and tth were adjusted as specified below. 

Some features of this model has already been discussed pre- 
viously H [U [T2I1 . The relevant solutions remain valid even 
for a 7^ 7. These solutions are the six homogeneous dis- 
tributions, the two cyclic defensive alliances, and three well 
mixed phases of two neutral species. For the cyclic defensive 
alliances the odd (or even) label species invade cyclically each 
other in the same wa y as it is described by the spatial Rock- 
Scissors-Paper game lIlSLfTilfTsll . The distinguished role (and 
also the name) of the cyclic defensive alliances come from the 
fact that the self-organizing spatio-temporal pattern provides 
a protection (stabihty) against external invaders liTi Im ITtIi . 

When X is increased for a = 7 = 1 the present system ex- 
hibits a first-order phase transition atX = Xc(Q! = 7= l) = 
0.05592(1) i^.lfX< Xc{a = 7 = 1) then one of the two 
cyclic defensive alliances will prevail the whole system after 
a domain growing process. Henceforth this final state will be 
denoted by referring to cyclic triplets. This model has 
three other defensive alliances composed from a neutral pair 



of species (e.g., and 3) because in their well-mixed phase the 
participants protect each other mutually against any external 
invaders 1^. The MC simulations have justified that one of 
these two-species defensive alliances will dominate the whole 
lattice after a domain growing process if X > Xc- This lat- 
ter final state is named in short as D (duet). Notice that only 
a portion of the species remains alive in this system for the 
uniform invasion rates. 

The internal symmetry of the two cycUc defensive alliances 
is conserved in the systems for the alliance-specific invasion 
rates li 1 lil . In the latter case four different invasion rates (a, 
(3, 7, and S) has been distinguished on the same food web 
plotted in Fig. [T] For X = 0, this type of parametrization has 
allowed us to study the cases where one of the cyclic defensive 
alliances is prefeiTed to the other. It turned out, for example, 
that the protection mechanism is enforced if the invasion rates 
are increased within a cyclic three-state alliance. This four- 
parameter model becomes equivalent to the present model for 
a — f3 and 7 = ^ in the absence of mixing. 

For the case of a = 1 and 7 = the food web has only one 
(six-species) cycle. This system was akeady investigated pre- 
viously by several authors llT2ilT8ll . In analogy to the spatial 
Rock-Scissors-Paper games the species alternates cyclically 
at each site and a self-organizing pattern is maintained by the 
moving invasion fronts for X = 0. 

For strong mixing the formation of well-mixed phases of 
the neutral species is expected. The three two-species asso- 
ciations are equivalent for a = 7 and the motion of inter- 
faces separating them is controlled by the curvature and ran- 
dom events |19]. This means that if two different domains 
are separated by a straight-line boundary then the average ve- 
locity of this interface is zero. However, if a > 7 then the 
well mixed phase of species and 3 can invade the teiTitory 
of the well-mixed phase of species 1 and 4, that can also in- 
vade the third association (consisting of species 2 and 5). In 
other words, the present model exemplifies a system where 
three associations play the spatial Rock-Scissors-Paper game. 
In the opposite case (a < 7) the direction of cyclic domi- 
nance is reversed. When visualizing the evolution of species 
distribution in this phase, one can recognize rotating spiral 
arms reported for maiiy other systems (for nice snapshots see 
the papers IH [T^ H^, IIH I22I1 and further references therein). 
This phase is denoted by T'^'{D) signaling the cyclic domi- 
nance of duets. Originally the recent research was planned 
to explore this phenomenon. It turned out, however, that the 
present model exhibits other self-organizing patterns as it will 
be detailed in the next section. 



III. THE RESULTS 

Without loosing generality we discuss separately the cases 
a < 7 (at 7 = 1) and 7 < a = 1. 
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A. The region a < 7 

First we study MC results obtained when varying X for 
a — 0.4 and 7 = 1. As previously discussed, the variation 
in the spatial distribution can be quantified by the above men- 
tioned pair correlation functions, namely, p„, pa, and p.y. In 
the numerical results plotted in Fig.|2]two arrows indicate the 
threshold values of the mixing {Xci(a) and Xc2{oi)) where 
phase transitions occur 

If X < Xci{a) then the finite system evolves into one 
of the T'-^ phases after a suitable relaxation (domain growth) 
time increasing with N. Within this phase the odd (or even) 
label species form a cyclic defensive alliance where the three 
mentioned species are present with the same average density 
(1/3) andp„ = = 0. 
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FIG. 2: The pair configuration probabilities pa (open squares), p.^ 
(closed squares) and p„ (open circles) as a function of X at fixed 
a = 0.4 and 7 = 1 values. Arrows show the positions of phase 
transitions. 

If X > Xc2{ct) then the three well mixed associations of 
the neutral pairs form a self-organizing pattern (T'~''{D)) re- 
sembling to the spatial Rock-Scissors-Paper game at a higher 
level. The typical extensions of domains and the width of 
boundary layers (separating two associations of neutral pairs) 
depend on X and a. The qualitative analysis indicates an in- 
crease in the typical domain size if a goes to 7 = 1 [providing 
X > Xci{a = 7=1) = Xc2{a = 7 = 1) = 0.05592(1)]. 
In fact, the driving force of the cyclic dominance is propor- 
tional to 7 — a. The numerical study of the impact of the 
vanishing cycUc dominance on the spatial distributions was 
akeady presented in a model combining the three-state Potts 
model and spatial Rock-Scissors-Paper game |23, 24]. In the 
light of the latter results it is expected that the typical domain 
size increases as l/|a — 7I when approaching the symmetric 
case (a = 7). 

The appearance of an intermediate region [Xd (a) < X < 
Xc2{ct)] in Fig. |2]was unexpected. The visualization of the 
evolution of species distribution (for a snapshot see Fig|3]l 
have indicated clearly that within this parameter region five 
types of domains (associations) can be distinguished. Namely, 



the two cyclic triplets (T) and also the three associations of 
neutral pairs (D) which form a self-organizing domain struc- 
ture. 
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FIG. 3: (Color online) Typical spatial distribution of species within 
a box of size 400 x 400 for a = 0.4, 7 = 1, and X = 0.02. 



It is worth emphasizing that a sufficiently large system 
size and long runs were necessary in the MC simulations to 
observe this intermediate region. More precisely, the self- 
organizing patterns has reached their final features (domain 
size, etc.) after a typical time of tth = 4 x 10^ MCSs if 
L — 4000. For the sake of comparison, the quantitative fea- 
tures of the T^{D) pattern can be well studied for L = 400 
after Uh ^ 4: x 10^ MCSs. 

Previous analyzes of similar systems have justified that the 
value of the critical point can also be determined by evaluating 
the average velocity of a straight invasion front separating two 
phases characterizing the final behavior below and above the 
critical point. The average velocity of this invasion front be- 
comes zero at the critical point. To clarify the behaviors in the 
intermediate region we have performed these numerical inves- 
tigations for different values of X. The results have clearly 
indicated that each D state can invade the territory of the T 
associations within the intermediate state. In other words, if 
an island of D (with a sufficiently large extension) is created 
via a nucleation mechanism within the territory of T (or even 
at the boundary of two T states) then this island grows perma- 
nently. 

In the present case, however, three equivalent D associ- 
ations exist which dominate cyclically each other as men- 
tioned above. Consequently, within the intermediate phase 
two growing D phases can collide and one of them will in- 
vade the other During the invasions the moving invasion 
front leave behind a slowly varying structure that differs from 
its final (well-mixed) distribution. Thus, the expanding D 
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FIG. 4: Variation of pair configuration probabilities [p-,{t) (solid) 
and p„ (t) (dashed line)] in the stationary state within the interme- 
diate region {a = 0.4, 7 = 1, and X = 0.021) for L = 4000. 
The values of pn are increased by a constant for easier comparison. 
Arrows show when D domains start to expand. 
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FIG. 5: Phase diagram of the model as a function of X and a for 7 = 
1. Ther^(s) region is characterized by the exclusive dominance of 
(0+2+4) or (1+3+5) cyclic alliance. The area T'^ (D) corresponds 
to the phase where three alliances of two-neutral species [(0 + 3), 
(1 + 4), and (2 + 5)] play spatial Rock-Scissors-Paper game. Within 
the shadowed region all the mentioned associations coexist and form 
a self-organizing pattern described in the text. 



associations become less stable against the invasion of the 
neighboring T associations in the vicinity of the moving in- 
vasion fronts. The visualization of the species distribution has 
demonstrated clearly that in many cases the newly invaded D 
territories were occupied by the neighboring T associations 
within a transient time. The alternative expansion of D and 
T domains can be traced well by monitoring the evolution of 
pair configuration probabilities. Evidently, the growth of D 
domains involves the increase of p„ while the extension of 
T domains increases the average value of p^. Consequently, 
one can observe opposite variations in the time-dependence 
of Pn{t) and Pj(t) as demonstrated in Fig.|4] Evidently, the 
"amplitude" of these variations vanish in the limit N ^ 00. 

Increasing X yields faster recovering (shorter transient 
time) and simultaneously makes the D territories more stable 
against the invasion of T domains. As a result, above a second 
threshold value (X > Xc2{ct)) the T associations cannot re- 
main alive and the whole system is prevailed by the previously 
described T*^ {D) phase. 

In order to determine the critical values of mixing (Xd (a) 
and Xc2{ct)), the MC simulations were repeated with increas- 
ing gradually the value of X for several values of a. The 
results are summarized in a phase diagram shown in Fig.|5] 

The intermediate region vanishes if a < ac — 0.170(5). 
More precisely, Xci{a) and Xc2{oi) go to zero simultane- 
ously if a tends to ac from above. For a < ac the T'^ [D) 
state occurs after a relaxation proportional to 1/X in the limit 
X ^ 0. In the absence of local mixing {X = 0), however, 
the well-mixed state of the neutral pairs cannot occur and the 
system develops into a state where the evolution of species 
distribution is governed by invasions of type 7 in a pattern 
exhibiting large domains of T associations. 



B. The region 7 < a 

Similarly to the previous section now we study the case of 
7 < a = 1 . Within this range of parameters the direction of 
cyclic invasions between the well mixed alliances of neutral 
species pair is reversed, that is, (0, 3) dominates (1, 4) domi- 
nates (2, 5) dominates (0, 3). As the direction of cyclic inva- 
sions [within the phase of T'^{D)] do not affect the the main 
features of pattern formation therefore we expect a phase di- 
agram similar to the case of a < 7. This expectation is sup- 
ported by Fig. |6] where the X dependences of pair configu- 
ration probabilities are plotted for 7 — 0.6. The qualitative 
similarity between Figs.|2]and|6]is striking. 
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FIG. 6: The pair configuration probabilities, such as pa (open 
squares), p^ (closed squares) and p„ (open circles) as a function of 
X at fixed 7 = 0.6 and a = 1 values. Arrows point to the positions 
of phase transitions. 
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Figure |6]represents a situation where both types (a and 7) 
of invasions play a relevant role in the pattern formation. In 
the opposite case (when either a or 7 tends to zero) a relevant 
difference occurs in the behaviors as plotted in Fig.|2]compar- 
ing pn changes as a function of a and 7 for a fixed value of 
X. 
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FIG. 7: The a and 7 dependence of p„ (open circles) and pid (closed 
squares) at a fixed value of mixing (X — 0.08). The arrow indicates 
the minimum of p„. 

Figure |7] shows that p„ has a local minimum at a small 
value of 7 which is missing in the case of a < 7. In fact, 
for 7 = the system develops into a state [denoted as S''^(s) 
(cyclic sextet)] where the six species invades cyclically each 
other. Within this spatio-temporal pattern the site-exchange 
process becomes rare and cannot affect significantly the spa- 
tial distribution of species. On the other hand, we can observe 
a smooth transition from the state T'^'{D) to S^{s) which is 
accompanied with the suppression of D domains (yielding a 
relevant decrease in p„) and with an increase of Pa and pid 
when decreasing the ratio -f/a. All these processes together 
result in a minimum in p„ that used to define a phase boundary 
between these phases. 
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FIG. 8: Phase diagram of the model as a function of X and 7 for 
a — 1. The notation of phases are the same as for Fig[5] S'^ (s) 
refers to a spatio-temporal pattern in which the evolution is domi- 
nantly governed by the cyclic invasions of the six species. 



The J — X phase diagram is shown in Fig. [8] where the 
dashed line shows the value of parameters where the mini- 
mum occurs in p„. According to our numerical investigations 
the value of Xci{'j) and Xc2{'^) coincide in this phase dia- 
gram within a range of a (0.35(2) < a < 0.45(2)) and both 
quantities vanish if a < 0.35(2). 



IV. SUMMARY 

We have studied a six-species ecological model on a square 
lattice where different types of associations are formed from 
a portion of species existing in the whole ecological model. 
The investigation of the present model was inspired by pre- 
vious results exemplifying several ways how the cyclic dom- 
inance can occur between the associations characterized by 
their composition and spatio-temporal pattern. In most of the 
previous studies the number of parameters was reduced by in- 
troducing many symmetries. Now we wished to explore some 
further phenomena that yield the formation and competition of 
alliances in more realistic biological systems when the sym- 
metries are reduced in the invasion rates. More precisely we 
have studied the cases characterized by two invasion rates, 
a and 7, in a way preserving the internal symmetries of the 
cyclic triplets. 

The present model exhibit a wide range of behaviors in the 
final stationary states as summarized in two phase diagrams 
(see Figs. |5] and [8]). For example, if a ^ 7 and the site ex- 
change mechanism is sufficiently strong then one observes a 
self-organizing spatio-temporal pattern in which the three al- 
liances of neutral pairs dominate cyclically each other. Al- 
though similar self-organizing patterns are reported in other 
systems |@], the present one seems to be the simplest lattice 
predator-prey model where the mechanism of cyclic domi- 
nance can take place at two different levels. In addition to 
this feature we have also revealed an unexpected phase where 
both the domains of cyclic three-species alliances and the neu- 
tral two-species alliances can coexist. The existence of this 
intermediate phase is closely related to the emergence of dif- 
ferent time- and length-scales within the self-organizing pat- 
terns. We think that further reduction of the symmetries in the 
species specific invasion rates can yield more complex behav- 
iors and other uncovered mechanisms supporting the coexis- 
tence of different alUances of species. 
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